library(marginaleffects)

# Start log file for script -----------------------------------------------
TeachingDemos::txtStart("./log files/figureA27_28_blackwhite_prep.txt")

# Load Data ---------------------------------------------------------
# individual data
load("./data/impdat1920.Race.WeekWeight.rdata")
load("./data/impdat20.Race.WeekWeight.rdata")

impdat20.WeekWeight <- fastDummies::dummy_cols(impdat20.WeekWeight, select_columns = "week")

# Equalizing weeks with Nationscape
impdat20.WeekWeight$week2 <- as.numeric(impdat20.WeekWeight$week) + 1
impdat20.WeekWeight$week2[which(impdat20.WeekWeight$week2 == 53)] <-  52
impdat20.WeekWeight$week2 <- as.character(impdat20.WeekWeight$week2)
impdat20.WeekWeight$week2[which(nchar(impdat20.WeekWeight$week2) == 1)] <- paste0("0", impdat20.WeekWeight$week2[which(nchar(impdat20.WeekWeight$week2) == 1)])
impdat20.WeekWeight <- fastDummies::dummy_cols(impdat20.WeekWeight, select_columns = "week2")

impdat20.WeekWeight$week2_post <- "Pre-GF"
for(i in 22:52){
  impdat20.WeekWeight$week2_post[which(impdat20.WeekWeight[, paste0("week2_", i)] == 1)] <- paste("Week", i)
}


vars <- c("year", "week", "date_mdy", "ideo_mod", "ideo_con", "age_cat4", "edu_cat4", "race5", "sex",
          "census_region", "broughtwebsite2", "WeekWeights", "t_wb")
impdat19_sub <- subset(impdat1920.WeekWeight, select = vars)
impdat20_sub <- subset(impdat20.WeekWeight, select = vars)

impdat_comb <- rbind(impdat19_sub, impdat20_sub)

# Race
impdat_comb$race_black <- ifelse(impdat_comb$race5 == "Black", 1, 0)
impdat_comb$race_black[which(is.na(impdat_comb$race5))] <- NA


# Race
impdat20.WeekWeight$race_black <- ifelse(impdat20.WeekWeight$race5 == "Black", "Black", "Not Black")
impdat20.WeekWeight$race_black[which(is.na(impdat20.WeekWeight$race5))] <- NA



# Models ------------------------------------------------------------------
m_week1920_bw <- lm(t_wb ~ 
                      as.factor(week)*year*race_black + ideo_mod + ideo_con +
                      age_cat4 + edu_cat4 + sex +
                      census_region + broughtwebsite2,
                    data = impdat_comb, weights = WeekWeights,
                    subset = date_mdy != "2020-05-25" & (race5 == "White" | race5 == "Black"))
summary(m_week1920_bw)


# full pre-post
m1_prepost_bw <- lm(t_wb ~ post_gf*race_black + 
                      age_cat4 + edu_cat4 + ideo7_sc + sex +
                      census_region + broughtwebsite2,
                    data = impdat20.WeekWeight, weights = WeekWeights,
                    subset = race5 == "White" | race5 == "Black")
summary(m1_prepost_bw)

# Weekly after
m1_pre_weeklypost_bw <- lm(t_wb ~ week2_post*race_black + 
                             age_cat4 + edu_cat4 + ideo7_sc + sex +
                             census_region + broughtwebsite2,
                           data = impdat20.WeekWeight, weights = WeekWeights,
                           subset = race5 == "White" | race5 == "Black")
summary(m1_pre_weeklypost_bw)


# Prediction Data Sets -----------------------------------------------------
weeks <- names(table(m_week1920_bw$model$`as.factor(week)`))

pdat_wht <- data.frame(week = weeks,
                       year = c(rep(2019, length(weeks)), rep(2020, length(weeks))),
                       ideo_mod = 0,
                       ideo_con = 0,
                       race_white = 1,
                       race_black = 0,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

pdat_blk <- data.frame(week = weeks,
                       year = c(rep(2019, length(weeks)), rep(2020, length(weeks))),
                       ideo_mod = 0,
                       ideo_con = 0,
                       race_white = 0,
                       race_black = 1,
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")


# Plot --------------------------------------------------------------------
preds_wht <- predict(m_week1920_bw, pdat_wht, se.fit = T)
preds_blk <- predict(m_week1920_bw, pdat_blk, se.fit = T)

p_dat <- data.frame(year = c(rep("2019", length(weeks)), rep("2020", length(weeks))),
                    week = weeks,
                    pred = c(preds_wht$fit,preds_blk$fit),
                    se = c(preds_wht$se.fit, preds_blk$se.fit),
                    race = c(rep("White", length(weeks)*2), rep("Black", length(weeks)*2)))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,54,107,160)] <- "2020-01-01" # from 00 weeks
p_dat$pred[which(p_dat$week == "2020-11-23" & p_dat$year == "2020")] <- NA
p_dat$se[which(p_dat$week == "2020-11-23" & p_dat$year == "2020")] <- NA
p_dat$pred[which(p_dat$week == "2020-11-30" & p_dat$year == "2020")] <- NA
p_dat$se[which(p_dat$week == "2020-11-30" & p_dat$year == "2020")] <- NA

ggplot(p_dat, aes(x = week, y = pred, group = year)) +
  geom_point(size = 2, aes(shape = year, color = year),
             position = position_dodge(width = .5)) +
  geom_linerange(aes(ymin = pred - 1.65*se,
                     ymax = pred + 1.65*se,
                     color = year),
                 position = position_dodge(width = .5),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred, color = year), se = FALSE) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred, color = year), se = FALSE) +
  theme_bw() +
  facet_grid(race~., scales = "free_y") +
  scale_color_manual(values = c("grey", "black"), name = "") +
  scale_shape_manual(values = c(1, 16), name = "") +
  labs(y = "Predicted D-Scores",
       x = "") +
  scale_y_continuous(labels = numform::ff_num(zero = 0, digits = 2)) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 20),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 20),
        strip.background = element_blank())

bw_1920 <- ggplot(p_dat, aes(x = week, y = pred, group = year)) +
  geom_point(size = 2, aes(shape = year, color = year),
             position = position_dodge(width = .5)) +
  geom_linerange(aes(ymin = pred - 1.39*se,
                     ymax = pred + 1.39*se,
                     color = year),
                 position = position_dodge(width = .5),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred, color = year), se = FALSE) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred, color = year), se = FALSE) +
  theme_bw() +
  facet_grid(race~.) +
  scale_color_manual(values = c("grey", "black"), name = "") +
  scale_shape_manual(values = c(1, 16), name = "") +
  labs(y = "Predicted White vs. Black Favorability",
       x = "") +
  scale_x_date(breaks = unique(p_dat$week)[c(1, 6, 10, 15, 19, 23, 28, 32, 37, 41, 45, 50)],
               labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
                          "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"),
               minor_breaks = NULL) + 
  scale_y_continuous(labels = numform::ff_num(zero = 0, digits = 2)) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 20),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 16),
        strip.background = element_blank())

p_dat_bw_therm <- p_dat

# Pre-Post Regressions ----------------------------------------------------
# pre-post full
d <- plot_cme(m1_prepost_bw, "post_gf", "race_black", draw = F)
prepost <- data.frame(race = c("Whites", "Blacks"),
                      coef = d$estimate,
                      se = d$std.error)

# pre weekly post
dd <- plot_cme(m1_pre_weeklypost_bw, "week2_post", "race_black", draw = F)
dd_wht <- subset(dd, race_black == "Not Black")
dd_blk <- subset(dd, race_black == "Black")
weekly_outs <- data.frame(race = c(rep("Blacks", nrow(dd_blk)),
                                   rep("Whites", nrow(dd_wht))),
                          week = c(gsub(" - Pre-GF", "", dd_blk$contrast),
                                   gsub(" - Pre-GF", "", dd_wht$contrast)),
                          coef = c(dd_blk$estimate, dd_wht$estimate),
                          se = c(dd_blk$std.error, dd_wht$std.error))

# Combine to plot
prepost <- subset(prepost, select = c("race", "coef", "se"))
prepost$week <- "Post GF, May 26, 2020-Dec 31, 2020"

coef_out <- rbind(prepost, weekly_outs)

coef_out$race <- factor(coef_out$race,
                        levels = c("Whites", "Blacks"))
coef_out$var <- as.factor(coef_out$week)
coef_out$var <- factor(coef_out$var,
                       levels = rev(levels(coef_out$var)))
coef_out$indx <- c(1, 1, rep(0, nrow(weekly_outs)))
coef_out$indx <- ifelse(coef_out$indx == 1, "Pre-Post", "Weekly Indicators")

coef_out_bw <- coef_out

# End log file for script -------------------------------------------------
TeachingDemos::txtStop()